# author:    Yannick Stiller, A. Duer, Robert A. Huber
# contact:   robert.huber@sbg.ac.at
# file name: analyses.R
# Context:   Provides analyses printed in Stiller et al World Trade Review

#Empty memory
rm(list = ls())

#Loads necessary packages
library(tidyverse)

#Load data
df <- read_rds("./created_data/df.rds")

#define useful functions for later
'%!in%' <- function(x,y)!('%in%'(x,y))
vcov_eff <- function(mod){sandwich::vcovCL(mod, cluster = df$dst)}

# Filter only countries for which competitiveness data is available
comp_available <- df %>%
  filter(!is.na(competitiveness)) %>%
  dplyr::select(country_short) %>%
  unique()

df <- df %>%
  filter(country_short %in% comp_available$country_short) %>%
  dplyr::select(country_short, region_name, region_code, trade_support, consequences_wages, consequences_jobs, 
                age, age_num, female, edu2, edu3, economic_left_right, competitiveness, employment, lgnic, GDPpc) %>%
  mutate(dst = sub("\\_u.*", "", region_code),
         dst = sub("\\_r.*", "", dst),
         consequences_wages_reg = factor(consequences_wages,
                                         levels = c("Does not make a difference", "Decrease", "Increase")),
         consequences_jobs_reg = factor(consequences_jobs,
                                        levels = c("Does not make a difference", "Job losses", "Job creation")))

# Save data as used in analysis
saveRDS(df, file = "./created_data/df_analysis.rds")

# Function to extract quantities of interest for texreg 
source("prepTex.R")

# Ordinal regression ####
#this section runs the main regression models

b_jobs_ord <- MASS::polr(data=df, consequences_jobs ~ edu2 + (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

b_wages_ord <- MASS::polr(data=df, consequences_wages ~ edu2 + (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

i_jobs_ord <- MASS::polr(data=df, consequences_jobs ~ edu2 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

i_wages_ord <- MASS::polr(data=df, consequences_wages ~ edu2 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

# Table output #
# this section exports the regression results for table 1 and 2.

texreg::texreg(l = list(prepTex(b_jobs_ord), prepTex(b_wages_ord)),
               file = "./tables/table1.tex",
               caption = "Education and the perceived consequences of trade",
               caption.above = T,
               label = "tab1",
               omit.coef = "country",
               custom.model.names = c("Jobs", "Wages"),
               custom.coef.names = c("Education (Tertiary)", "Age (41-65)",
                                     "Age (66+)", "Logged Regional GNIpc", "Employment (Employed)",
                                     "Subnational Trade Competitiveness", "Economic Left-Right", "Gender (Female)"),
               custom.note = "\\parbox{0.75\\linewidth}{%stars. Entries are unstandardized coefficients from a logistic regression. Standard errors (in brackets) are clustered on a regional level. Country-fixed effects omitted from the table.}",
               single.row = T,
               stars = c(0.1, 0.05, 0.01),
               digits = 2,
               leading.zero = T, 
               reorder.coef = c(1, 5, 2:3, 4, 6:8),
               dcolumn = T,
               use.packages = F,
               center = T,
               float.pos = "htb")

texreg::texreg(l = list(prepTex(i_jobs_ord), prepTex(i_wages_ord)),
               file = "./tables/table2.tex",
               caption = "Education, economic context, and perceptions of the consequences of trade",
               caption.above = T,
               label = "tab2",
               omit.coef = "country",
               custom.model.names = c("Jobs", "Wages"), 
               custom.coef.names = c("Education (Tertiary)", "Age (41-65)",
                                     "Age (66+)", "Logged Regional GNIpc", "Employment (Employed)",
                                     "Subnational Trade Competitiveness", "Economic Left-Right", "Gender (Female)",
                                     "Edu. x Age (41-65)", "Edu. x Age (65+)",
                                     "Edu. x Dst. GNIpc", "Edu. x Employment", "Edu. x Comp."),
               custom.note = "\\parbox{0.75\\linewidth}{%stars. Entries are unstandardized coefficients from a logistic regression. Standard errors (in brackets) are clustered on a regional level. Country-fixed effects omitted from the table.}",
               single.row = T,
               stars = c(0.1, 0.05, 0.01),
               digits = 2,
               leading.zero = T, 
               reorder.coef = c(1, 5, 12, 2:3, 9:10, 4, 11, 6, 13, 7:8),
               dcolumn = T,
               use.packages = F,
               center = T,
               float.pos = "htb")

# Predictions #
#this section exports the central figures based on the aforementioned models for the appendix and maintext.

# .employment status ####

df_text <- data.frame(employment = rep(c(1.65, 1.95)),
                      pred = c(0.625, 0.58),
                      dim = factor(c("Trade induces job creation", "Trade induces job creation")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:employment", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu2:employment", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu2", "employment", "dim")) %>%
  ggplot(., aes(x=employment, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nEmployment Status", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA6.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(employment = rep(c(1.65, 1.95)),
                      pred = c(0.470, .425),
                      dim = factor(c("Trade increases wages", "Trade increases wages")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:employment", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu2:employment", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu2", "employment", "dim")) %>%
  ggplot(., aes(x=employment, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nEmployment Status", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA7.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(employment = rep(c(1.75, 1.6), 2),
                      pred = c(0.624, 0.56, 0.465, 0.445),
                      dim = factor(rep(c("Trade induces job creation", "Trade increases wages"), each = 2)),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

rbind(dplyr::left_join(as.data.frame(effects::effect("edu2:employment", i_jobs_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    employment = c("Unemployed", "Employed")),
                                                     x.var="employment", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, employment, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                                    ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                             levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                       as.data.frame(effects::effect("edu2:employment", i_jobs_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    employment = c("Unemployed", "Employed")),
                                                     x.var="employment", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, employment, contains("se.prob")) %>% 
                         gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                                    ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                             levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                       by = c("edu2", "employment", "dim")) %>%
        mutate(dv = "Jobs"),
      dplyr::left_join(as.data.frame(effects::effect("edu2:employment", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    employment = c("Unemployed", "Employed")),
                                                     x.var="employment", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, employment, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       as.data.frame(effects::effect("edu2:employment", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    employment = c("Unemployed", "Employed")),
                                                     x.var="employment", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, employment, contains("se.prob")) %>% 
                         gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       by = c("edu2", "employment", "dim")) %>%
        mutate(dv = "Wages")) %>%
  filter(dim %in% c("Trade induces job creation", "Trade increases wages")) %>%
  ggplot(., aes(x=employment, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nEmployment Status", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figure1.pdf",
       width = 20, height = 9, units = "cm")

# .Age ####

df_text <- data.frame(age = rep(c(2.5, 2.25)),
                      pred = c(0.675, 0.575),
                      dim = factor(c("Trade induces job creation", "Trade induces job creation")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:age", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu2:age", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu2", "age", "dim")) %>%
  ggplot(., aes(x=age, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge group", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA8.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(age = rep(c(2.5, 2.55)),
                      pred = c(0.5, 0.42),
                      dim = factor(c("Trade increases wages", "Trade increases wages")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:age", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu2:age", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu2", "age", "dim")) %>%
  ggplot(., aes(x=age, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge group", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA9.pdf",
       width = 20, height = 9, units = "cm")


df_text <- data.frame(age = rep(c(2.75, 2.5), 2),
                      pred = c(0.675, 0.585, 0.5, 0.43),
                      dim = factor(rep(c("Trade induces job creation", "Trade increases wages"), each = 2)),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

rbind(dplyr::left_join(as.data.frame(effects::effect("edu2:age", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu2:age", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu2", "age", "dim")) %>%
        mutate(dv = "Jobs"),
      dplyr::left_join(as.data.frame(effects::effect("edu2:age", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    age = c("young", "middle", "old")),
                                                     x.var="age", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, age, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       as.data.frame(effects::effect("edu2:age", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    age = c("young", "middle", "old")),
                                                     x.var="age", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, age, contains("se.prob")) %>% 
                         gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       by = c("edu2", "age", "dim")) %>%
        mutate(dv = "Wages")) %>%
  filter(dim %in% c("Trade induces job creation", "Trade increases wages")) %>%
  ggplot(., aes(x=age, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge group", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figure2.pdf",
       width = 20, height = 9, units = "cm")

# .level of development ####

df_text <- data.frame(lgnic = rep(c(3.25, 2.125)),
                      pred = c(0.64, 0.525),
                      dim = factor(c("Trade induces job creation", "Trade induces job creation")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))


dplyr::left_join(as.data.frame(effects::effect("edu2:lgnic", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu2:lgnic", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu2", "lgnic", "dim")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA10.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(lgnic = rep(c(3.25, 2.25)),
                      pred = c(.47, .56),
                      dim = factor(c("Trade increases wages", "Trade increases wages")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:lgnic", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu2:lgnic", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu2", "lgnic", "dim")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA11.pdf",
       width = 20, height = 9, units = "cm")


df_text <- data.frame(lgnic = rep(c(2.75, 2.45), 2),
                      pred = c(0.64, 0.51, 0.49, 0.39),
                      dim = factor(rep(c("Trade induces job creation", "Trade increases wages"), each = 2)),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

rbind(dplyr::left_join(as.data.frame(effects::effect("edu2:lgnic", i_jobs_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                                to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                                length.out = 11)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                                    ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                             levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                       as.data.frame(effects::effect("edu2:lgnic", i_jobs_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                                to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                                length.out = 11)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, lgnic, contains("se.prob")) %>% 
                         gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                                    ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                             levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                       by = c("edu2", "lgnic", "dim")) %>%
        mutate(dv = "Jobs"),
      dplyr::left_join(as.data.frame(effects::effect("edu2:lgnic", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                                to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                                length.out = 11)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       as.data.frame(effects::effect("edu2:lgnic", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                                to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                                length.out = 11)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, lgnic, contains("se.prob")) %>% 
                         gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       by = c("edu2", "lgnic", "dim")) %>%
        mutate(dv = "Wages")) %>%
  filter(dim %in% c("Trade induces job creation", "Trade increases wages")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figure3.pdf",
       width = 20, height = 9, units = "cm")

# .RCA ####

df_text <- data.frame(competitiveness = rep(c(1.25, 1.125)),
                      pred = c(0.7, 0.59),
                      dim = factor(c("Trade induces job creation", "Trade induces job creation")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:competitiveness", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu2:competitiveness", i_jobs_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, competitiveness, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu2", "competitiveness", "dim")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA12.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(competitiveness = rep(c(1.25, .95)),
                      pred = c(.545, .42),
                      dim = factor(c("Trade increases wages", "Trade increases wages")),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

dplyr::left_join(as.data.frame(effects::effect("edu2:competitiveness", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                             to = 2,
                                                                             length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu2:competitiveness", i_wages_ord, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                             to = 2,
                                                                             length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, competitiveness, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu2", "competitiveness", "dim")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA13.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(competitiveness = rep(c(0.65, 1.5), 2),
                      pred = c(0.66, 0.585, 0.51, 0.415),
                      dim = factor(rep(c("Trade induces job creation", "Trade increases wages"), each = 2)),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Tertiary", "Non-Tertiary"),2))

rbind(dplyr::left_join(as.data.frame(effects::effect("edu2:competitiveness", i_jobs_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    competitiveness = seq(from = -2,
                                                                                   to = 2,
                                                                                   length.out = 9)),
                                                     x.var="competitiveness", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                                    ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                             levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                       as.data.frame(effects::effect("edu2:competitiveness", i_jobs_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    competitiveness = seq(from = -2,
                                                                                   to = 2,
                                                                                   length.out = 9)),
                                                     x.var="competitiveness", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, competitiveness, contains("se.prob")) %>% 
                         gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                                    ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                             levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                       by = c("edu2", "competitiveness", "dim")) %>%
        mutate(dv = "Jobs"),
      dplyr::left_join(as.data.frame(effects::effect("edu2:competitiveness", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    competitiveness = seq(from = -2,
                                                                                   to = 2,
                                                                                   length.out = 9)),
                                                     x.var="competitiveness", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       as.data.frame(effects::effect("edu2:competitiveness", i_wages_ord, vcov = vcov_eff,
                                                     xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                                    competitiveness = seq(from = -2,
                                                                                   to = 2,
                                                                                   length.out = 9)),
                                                     x.var="competitiveness", confidence.level = 0.9)) %>% 
                         dplyr::select(edu2, competitiveness, contains("se.prob")) %>% 
                         gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       by = c("edu2", "competitiveness", "dim")) %>%
        mutate(dv = "Wages")) %>%
  filter(dim %in% c("Trade induces job creation", "Trade increases wages")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figure4.pdf",
       width = 20, height = 9, units = "cm")


# Robustness: Trade Support ####
#This section runs the models for section 5.3 and the respective tables

b_trade_ord <- MASS::polr(data=df, trade_support ~ edu2 + (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)
b_trade_ord_dc <- MASS::polr(data=df, trade_support ~ consequences_jobs_reg + consequences_wages_reg + edu2 + (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

texreg::texreg(l = list(prepTex(b_trade_ord), prepTex(b_trade_ord_dc)),
               file = "./tables/tableA3.tex",
               caption = "Education, Distributional Consequences of Trade, and Trade Support",
               caption.above = T,
               label = "tabA3",
               omit.coef = "country",
               custom.model.names = c("Trade Support", "Trade Support"),
               custom.coef.names = c("Education (Tertiary)", "Age (41-65)",
                                     "Age (66+)", "Logged Regional GNIpc", "Employment (Employed)",
                                     "Subnational Trade Competitiveness", "Economic Left-Right", "Gender (Female)",
                                     "Trade induces job losses", "Trade induces job creation",
                                     "Trade decreases wages", "Trade increases wages"),
               custom.note = "\\parbox{0.85\\linewidth}{%stars. Entries are unstandardized coefficients from a logistic regression. Standard errors (in brackets) are clustered on a regional level. Country-fixed effects omitted from the table.}",
               single.row = T,
               stars = c(0.1, 0.05, 0.01),
               digits = 2,
               leading.zero = T,
               reorder.coef = c(1, 9:12, 2:3, 5, 4, 6:8),
               dcolumn = T,
               use.packages = F,
               center = T,
               float.pos = "htb")

# Robustness: Education 3 #### 
#This section provide the analyses for section A7.1 in the Online Appendix

b_jobs_ord3 <- MASS::polr(data=df, consequences_jobs ~ edu3 + (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

b_wages_ord3 <- MASS::polr(data=df, consequences_wages ~ edu3 + (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

i_jobs_ord3 <- MASS::polr(data=df, consequences_jobs ~ edu3 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

i_wages_ord3 <- MASS::polr(data=df, consequences_wages ~ edu3 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

# .employment status ####

df_text <- data.frame(employment = c(1.53, 2, 1.75),
                      pred = c(0.54, 0.58, 0.625),
                      dim = factor(rep("Trade induces job creation", 3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu3:employment", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, employment, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu3:employment", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, employment, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu3", "employment", "dim")) %>%
  ggplot(., aes(x=employment, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .5)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .5), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nEmployment Status", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA16.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(employment = c(1.525, 2.3, 1.525),
                      pred = c(0.447, 0.444, 0.467),
                      dim = factor(rep("Trade increases wages", 3)),
                      edu3 = c("Primary", "Sec.", "Tertiary"),
                      lab = c("Primary", "Sec.", "Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu3:employment", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, employment, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu3:employment", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed")),
                                               x.var="employment", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, employment, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu3", "employment", "dim")) %>%
  mutate(edu3 = factor(edu3,
                       levels = c("Primary", "Secondary", "Tertiary"),
                       labels = c("Primary", "Sec.", "Tertiary"))) %>%
  ggplot(., aes(x=employment, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .25)) + 
  geom_text(data = df_text, mapping = aes(label = lab), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nEmployment Status", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA17.pdf",
       width = 20, height = 9, units = "cm")

# .Age ####

df_text <- data.frame(age = rep(c(2.5, 2, 2.5)),
                      pred = c(0.525, 0.59, 0.675),
                      dim = factor(rep(c("Trade induces job creation"),3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu3:age", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, age, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu3:age", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, age, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu3", "age", "dim")) %>%
  ggplot(., aes(x=age, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .5)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .5), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge group", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA18.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(age = c(2.5, 1.5, 2.5),
                      pred = c(0.425, 0.48, 0.51),
                      dim = factor(rep(c("Trade increases wages"),3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu3:age", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, age, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu3:age", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              age = c("young", "middle", "old")),
                                               x.var="age", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, age, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu3", "age", "dim")) %>%
  ggplot(., aes(x=age, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .5)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .5), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge group", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA19.pdf",
       width = 20, height = 9, units = "cm")

# .level of development ####

df_text <- data.frame(lgnic = c( 3.1, 2.5, 3.25),
                      pred = c(0.415, 0.59, .65),
                      dim = factor(rep("Trade induces job creation", 3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary")) %>%
  filter(edu3 != "Secondary")

dplyr::left_join(as.data.frame(effects::effect("edu3:lgnic", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu3:lgnic", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu3", "lgnic", "dim")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA20.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(lgnic = c( 3.25, 1.75, 3.25),
                      pred = c(0.425, 0.605, .47),
                      dim = factor(rep("Trade increases wages", 3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary")) %>%
  filter(edu3 != "Primary")


dplyr::left_join(as.data.frame(effects::effect("edu3:lgnic", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu3:lgnic", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu3", "lgnic", "dim")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA21.pdf",
       width = 20, height = 9, units = "cm")

# .RCA ####

df_text <- data.frame(competitiveness = c( 0, rep(1.25, 2)),
                      pred = c(0.5, 0.6, .7),
                      dim = factor(rep("Trade induces job creation", 3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu3:competitiveness", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu3:competitiveness", i_jobs_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, competitiveness, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu3", "competitiveness", "dim")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA22.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(competitiveness = c( 0, -1, 1.25),
                      pred = c(0.425, 0.49, .55),
                      dim = factor(rep("Trade increases wages", 3)),
                      edu3 = c("Primary", "Secondary", "Tertiary"),
                      lab = c("Primary", "Secondary", "Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu3:competitiveness", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu3:competitiveness", i_wages_ord3, vcov = vcov_eff,
                                               xlevels = list(edu3 = c("Primary", "Secondary", "Tertiary"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu3, competitiveness, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu3", "competitiveness", "dim")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = edu3, shape = edu3)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA23.pdf",
       width = 20, height = 9, units = "cm")

# Robustness: Age as continuous #### 
#This section provide the analyses for Appendix section 7.2

i_jobs_ord_age_con <- MASS::polr(data=df, consequences_jobs ~ edu2 * (age_num + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)
i_wages_ord_age_con <- MASS::polr(data=df, consequences_wages ~ edu2 * (age_num + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

df_text <- data.frame(age_num = rep(c(35,35)),
                      pred = c(0.62, 0.54),
                      dim = factor(rep("Trade induces job creation", 2)),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = c("Tertiary", "Non-Tertiary"))


dplyr::left_join(as.data.frame(effects::effect("edu2:age_num", i_jobs_ord_age_con, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age_num = seq(from = 20,
                                                                          to = 80,
                                                                          length.out = 11)),
                                               x.var="age_num", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age_num, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("edu2:age_num", i_jobs_ord_age_con, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age_num = seq(from = 20,
                                                                          to = 80,
                                                                          length.out = 11)),
                                               x.var="age_num", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age_num, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 by = c("edu2", "age_num", "dim")) %>%
  ggplot(., aes(x=age_num, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = 2)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = 2), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA24.pdf",
       width = 20, height = 9, units = "cm")

df_text <- data.frame(age_num = rep(c(55,55)),
                      pred = c(0.52, 0.415),
                      dim = factor(rep("Trade increases wages", 2)),
                      edu2 = c("Tertiary", "Non-Tertiary"),
                      lab = c("Tertiary", "Non-Tertiary"))

dplyr::left_join(as.data.frame(effects::effect("edu2:age_num", i_wages_ord_age_con, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age_num = seq(from = 20,
                                                                          to = 80,
                                                                          length.out = 11)),
                                               x.var="age_num", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age_num, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 as.data.frame(effects::effect("edu2:age_num", i_wages_ord_age_con, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              age_num = seq(from = 20,
                                                                          to = 80,
                                                                          length.out = 11)),
                                               x.var="age_num", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, age_num, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                       levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                 by = c("edu2", "age_num", "dim")) %>%
  ggplot(., aes(x=age_num, y = pred, colour = edu2, shape = edu2)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = 2)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = 2), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nAge", y = "Predicted probability\n") +
  scale_colour_discrete("Education:") + 
  scale_shape_discrete("Education:") + 
  theme(legend.position = "bottom") + 
  facet_wrap(dim~., scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA25.pdf",
       width = 20, height = 9, units = "cm")

# Robustness: Threeway interaction #### 
#This section provides the analyses for section A7.3 of the online appendix

i_jobs_ord_3w_rca <- MASS::polr(data=df, consequences_jobs ~ edu2 * employment * competitiveness + age + lgnic + economic_left_right + female + country_short, Hess = T)
i_wages_ord_3w_rca <- MASS::polr(data=df, consequences_wages ~ edu2 * employment * competitiveness + age + lgnic + economic_left_right + female + country_short, Hess = T)

i_jobs_ord_3w_gni <- MASS::polr(data=df, consequences_jobs ~ edu2 * employment * lgnic + age +  competitiveness + economic_left_right + female + country_short, Hess = T)
i_wages_ord_3w_gni <- MASS::polr(data=df, consequences_wages ~ edu2 * employment * lgnic + age + competitiveness + economic_left_right + female + country_short, Hess = T)



# . level of development ####

df_text <- data.frame(lgnic = c( 3, 2.5),
                      employment = rep(c("Employed", "Unemployed"),each = 2),
                      pred = rep(c(0.65, 0.525), each = 2),
                      dim = factor(c("Trade induces\njob creation", "Trade induces\njob creation")),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Employed", "Unemployed"),each = 2)) %>%
  filter(edu2 != "Non-Tertiary")

dplyr::left_join(as.data.frame(effects::effect("edu2:employment:lgnic", i_jobs_ord_3w_gni, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces\njob losses",
                                              ifelse(grepl("crea", dim), "Trade induces\njob creation", "Trade does not\nmake a difference")),
                                       levels = c("Trade induces\njob losses", "Trade does not\nmake a difference", "Trade induces\njob creation"))),
                 as.data.frame(effects::effect("edu2:employment:lgnic", i_jobs_ord_3w_gni, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces\njob losses",
                                              ifelse(grepl("crea", dim), "Trade induces\njob creation", "Trade does not\nmake a difference")),
                                       levels = c("Trade induces\njob losses", "Trade does not\nmake a difference", "Trade induces\njob creation"))),
                 by = c("edu2", "employment", "lgnic", "dim")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = employment, shape = employment)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Employment:") + 
  scale_shape_discrete("Employment:") + 
  theme(legend.position = "bottom") + 
  facet_grid(dim~edu2, scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA26.pdf",
       width = 20, height = 15, units = "cm")


df_text <- data.frame(lgnic = c(3, 2.0, 2.5, 2),
                      employment = rep(c("Employed", "Unemployed"),each = 2),
                      pred = rep(c(0.5, 0.4), each = 2),
                      dim = factor(c("Trade increases\nwages", "Trade increases\nwages")),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Employed", "Unemployed"),each = 2)) %>%
  filter(edu2 != "Non-Tertiary")

dplyr::left_join(as.data.frame(effects::effect("edu2:employment:lgnic", i_wages_ord_3w_gni, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases\nwages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases\nwages", "Trade does not\nmake a difference")),
                                       levels = c("Trade decreases\nwages", "Trade does not\nmake a difference", "Trade increases\nwages"))),
                 as.data.frame(effects::effect("edu2:employment:lgnic", i_wages_ord_3w_gni, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              lgnic = seq(from = quantile(df$lgnic, 0.1, na.rm = T),
                                                                          to = quantile(df$lgnic, 0.9, na.rm = T),
                                                                          length.out = 11)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment,lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases\nwages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases\nwages", "Trade does not\nmake a difference")),
                                       levels = c("Trade decreases\nwages", "Trade does not\nmake a difference", "Trade increases\nwages"))),
                 by = c("edu2", "employment", "lgnic", "dim")) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = employment, shape = employment)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Employment:") + 
  scale_shape_discrete("Employment:") + 
  theme(legend.position = "bottom") + 
  facet_grid(dim~edu2, scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA27.pdf",
       width = 20, height = 15, units = "cm")

# .RCA ####

df_text <- data.frame(competitiveness = rep(c(1.275, 1.125),2),
                      employment = rep(c("Employed", "Unemployed"),each = 2),
                      pred = rep(c(0.74, 0.575), each = 2),
                      dim = factor(c("Trade induces\njob creation", "Trade induces\njob creation")),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Employed", "Unemployed"),each = 2)) %>%
  filter(edu2 != "Non-Tertiary")

dplyr::left_join(as.data.frame(effects::effect("edu2:employment:competitiveness", i_jobs_ord_3w_rca, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces\njob losses",
                                              ifelse(grepl("crea", dim), "Trade induces\njob creation", "Trade does not\nmake a difference")),
                                       levels = c("Trade induces\njob losses", "Trade does not\nmake a difference", "Trade induces\njob creation"))),
                 as.data.frame(effects::effect("edu2:employment:competitiveness", i_jobs_ord_3w_rca, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment,competitiveness, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces\njob losses",
                                              ifelse(grepl("crea", dim), "Trade induces\njob creation", "Trade does not\nmake a difference")),
                                       levels = c("Trade induces\njob losses", "Trade does not\nmake a difference", "Trade induces\njob creation"))),
                 by = c("edu2", "employment", "competitiveness", "dim")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = employment, shape = employment)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Employment:") + 
  scale_shape_discrete("Employment:") + 
  theme(legend.position = "bottom") + 
  facet_grid(dim~edu2, scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA28.pdf",
       width = 20, height = 15, units = "cm")


df_text <- data.frame(competitiveness = rep(c(1.275, 1.275),each = 2),
                      employment = rep(c("Employed", "Unemployed"),each = 2),
                      pred = rep(c(0.6, 0.415), each = 2),
                      dim = factor(c("Trade increases\nwages", "Trade increases\nwages")),
                      edu2 = rep(c("Tertiary", "Non-Tertiary"),2),
                      lab = rep(c("Employed", "Unemployed"),each = 2)) %>%
  filter(edu2 != "Non-Tertiary")

dplyr::left_join(as.data.frame(effects::effect("edu2:employment:competitiveness", i_wages_ord_3w_rca, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment, competitiveness, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases\nwages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases\nwages", "Trade does not\nmake a difference")),
                                       levels = c("Trade decreases\nwages", "Trade does not\nmake a difference", "Trade increases\nwages"))),
                 as.data.frame(effects::effect("edu2:employment:competitiveness", i_wages_ord_3w_rca, vcov = vcov_eff,
                                               xlevels = list(edu2 = c("Non-Tertiary", "Tertiary"),
                                                              employment = c("Unemployed", "Employed"),
                                                              competitiveness = seq(from = -2,
                                                                                     to = 2,
                                                                                     length.out = 9)),
                                               x.var="competitiveness", confidence.level = 0.9)) %>% 
                   dplyr::select(edu2, employment,competitiveness, contains("se.prob")) %>% 
                   gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases\nwages",
                                              ifelse(grepl("ecrease", dim), "Trade decreases\nwages", "Trade does not\nmake a difference")),
                                       levels = c("Trade decreases\nwages", "Trade does not\nmake a difference", "Trade increases\nwages"))),
                 by = c("edu2", "employment", "competitiveness", "dim")) %>%
  ggplot(., aes(x=competitiveness, y = pred, colour = employment, shape = employment)) + 
  geom_pointrange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  geom_text(data = df_text, mapping = aes(label = lab), position = position_dodge(width = .1), show.legend = FALSE) + 
  theme_minimal() + 
  labs(x = "\nSubnational Trade Competitiveness", y = "Predicted probability\n") +
  scale_colour_discrete("Employment:") + 
  scale_shape_discrete("Employment:") + 
  theme(legend.position = "bottom") + 
  facet_grid(dim~edu2, scales = "free_y") + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA29.pdf",
       width = 20, height = 15, units = "cm")

# Split by development ####
#this section provides analyses reported in Figures A14 and A15

source("prepTex_split.R")

vcov_effH <- function(mod){sandwich::vcovCL(mod, cluster = df_H$dst)}
vcov_effL <- function(mod){sandwich::vcovCL(mod, cluster = df_L$dst)}

df_L <- df %>% filter(GDPpc <= mean(df$GDPpc, na.rm=T))
df_H <- df %>% filter(GDPpc > mean(df$GDPpc, na.rm=T))


i_jobs_ordH <- MASS::polr(data=df_H, consequences_jobs ~ edu2 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)
i_jobs_ordL <- MASS::polr(data=df_L, consequences_jobs ~ edu2 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

i_wages_ordH <- MASS::polr(data=df_H, consequences_wages ~ edu2 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)
i_wages_ordL <- MASS::polr(data=df_L, consequences_wages ~ edu2 * (age + lgnic + employment + competitiveness) + economic_left_right + female + country_short, Hess = T)

rbind(dplyr::left_join(as.data.frame(effects::effect("lgnic", i_jobs_ordH, vcov = vcov_effH,
                                               xlevels = list(lgnic = seq(from = quantile(df_H$lgnic, 0.0, na.rm = T),
                                                                          to = quantile(df_H$lgnic, 1, na.rm = T),
                                                                          length.out = 7)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("lgnic", i_jobs_ordH, vcov = vcov_effH,
                                               xlevels = list(
                                                              lgnic = seq(from = quantile(df_H$lgnic, 0, na.rm = T),
                                                                          to = quantile(df_H$lgnic, 1, na.rm = T),
                                                                          length.out = 7)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation")),
                          dev = "High"),
                 by = c("lgnic", "dim")),
      dplyr::left_join(as.data.frame(effects::effect("lgnic", i_jobs_ordL, vcov = vcov_effL,
                                               xlevels = list(
                                                 lgnic = seq(from = quantile(df_L$lgnic, 0, na.rm = T),
                                                             to = quantile(df_L$lgnic, 1, na.rm = T),
                                                             length.out = 7)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                   gather(prob.Job.losses:prob.Job.creation, key = "dim", value = "pred") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation"))),
                 as.data.frame(effects::effect("lgnic", i_jobs_ordL, vcov = vcov_effL,
                                               xlevels = list(
                                                 lgnic = seq(from = quantile(df_L$lgnic, 0, na.rm = T),
                                                             to = quantile(df_L$lgnic, 1, na.rm = T),
                                                             length.out = 7)),
                                               x.var="lgnic", confidence.level = 0.9)) %>% 
                   dplyr::select(lgnic, contains("se.prob")) %>% 
                   gather(se.prob.Job.losses:se.prob.Job.creation, key = "dim", value = "se") %>%
                   mutate(dim = factor(ifelse(grepl("loss", dim), "Trade induces job losses",
                                              ifelse(grepl("crea", dim), "Trade induces job creation", "Trade does not make a difference")),
                                       levels = c("Trade induces job losses", "Trade does not make a difference", "Trade induces job creation")),
                          dev = "Low"),
                 by = c("lgnic", "dim"))) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = dev)) + 
  geom_linerange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_point(size = 1, position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Development (country):") + 
  theme(legend.position = "bottom") + 
  facet_wrap(.~dim) + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA14.pdf",
       width = 20, height = 9, units = "cm")

rbind(dplyr::left_join(as.data.frame(effects::effect("lgnic", i_wages_ordH, vcov = vcov_effH,
                                                     xlevels = list(
                                                       lgnic = seq(from = quantile(df_H$lgnic, 0, na.rm = T),
                                                                   to = quantile(df_H$lgnic, 1, na.rm = T),
                                                                   length.out = 7)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       as.data.frame(effects::effect("lgnic", i_wages_ordH, vcov = vcov_effH,
                                                     xlevels = list(
                                                       lgnic = seq(from = quantile(df_H$lgnic, 0, na.rm = T),
                                                                   to = quantile(df_H$lgnic, 1, na.rm = T),
                                                                   length.out = 7)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(lgnic, contains("se.prob")) %>% 
                         gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages")),
                                dev = "High"),
                       by = c("lgnic", "dim")),
      dplyr::left_join(as.data.frame(effects::effect("lgnic", i_wages_ordL, vcov = vcov_effL,
                                                     xlevels = list(
                                                       lgnic = seq(from = quantile(df_L$lgnic, 0, na.rm = T),
                                                                   to = quantile(df_L$lgnic, 1, na.rm = T),
                                                                   length.out = 7)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(lgnic, contains("prob"), -contains("se."), -contains("L."), -contains("U.")) %>%
                         gather(prob.Decrease:prob.Increase, key = "dim", value = "pred") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages"))),
                       as.data.frame(effects::effect("lgnic", i_wages_ordL, vcov = vcov_effL,
                                                     xlevels = list(
                                                       lgnic = seq(from = quantile(df_L$lgnic, 0, na.rm = T),
                                                                   to = quantile(df_L$lgnic, 1, na.rm = T),
                                                                   length.out = 7)),
                                                     x.var="lgnic", confidence.level = 0.9)) %>% 
                         dplyr::select(lgnic, contains("se.prob")) %>% 
                         gather(se.prob.Decrease:se.prob.Increase, key = "dim", value = "se") %>%
                         mutate(dim = factor(ifelse(grepl("ncrease", dim), "Trade increases wages",
                                                    ifelse(grepl("ecrease", dim), "Trade decreases wages", "Trade does not make a difference")),
                                             levels = c("Trade decreases wages", "Trade does not make a difference", "Trade increases wages")),
                                dev = "Low"),
                       by = c("lgnic", "dim"))) %>%
  ggplot(., aes(x=lgnic, y = pred, colour = dev)) + 
  geom_linerange(aes(ymin = pred - 1.645*se, ymax = pred + 1.645*se), position = position_dodge(width = .1)) + 
  geom_point(size = 1, position = position_dodge(width = .1)) + 
  geom_line(alpha = .2, position = position_dodge(width = .1)) + 
  theme_minimal() + 
  labs(x = "\nLogged Regional GNI per Capita", y = "Predicted probability\n") +
  scale_colour_discrete("Development (country):") + 
  theme(legend.position = "bottom") + 
  facet_wrap(.~dim) + 
  theme(strip.placement = "outside",
        panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  NULL

ggsave("./figures/figureA15.pdf",
       width = 20, height = 9, units = "cm")
